In this paper, we use the NOURISHING framework and policy database to assess the implementation of healthy eating policies across Europe.
There is limited knowledge of implemented nutrition policies in Europe. The aim of this report is to showcase index construction by providing an overview of implemented policies in different European countries, and identify any similarities or clusters of policies between European countries. A simple regression displaying the index against obesity rates in European countries is included.
The findings indicate an overall low implementation of nutrition policies across European countries using the NOURISHING framework. A country may have multiple policies but an holistic approach to the policy landscape is missing.
For future index construction, the low policy coverage may pose a problem for an index using the benchmarks. This can affect aggregation, hence a scoreboard may be a better approach than an overall index. More investigation is necessary: correlation between aggregation levels, principal components analyses and sensitivity analyses.
The data is manipulated in R prior to analyses. The following packages are utilised:
Data is open source from http://policydatabase.wcrf.org.
policydb <- read_csv("C:\\Users\\JRMA\\OneDrive - Folkehelseinstituttet\\Dokumenter\\CO-CREATE\\policy-export10-Dec-2021.csv", show_col_types = FALSE)
colnames(policydb) <- c("DB_type", "PolicyArea", "SubPArea", "PolicyAction", "Country", "Topics", "BenchmarkID", "PolicyAreaID", "PolicyDim")
policydb <- policydb[-c(1,6)]
country_hierch2 <- policydb %>% count(Country, PolicyAreaID, BenchmarkID)
truef_spa <- country_hierch2 %>% mutate(
PolicyActionImplemented = if_else(n > 0, TRUE, FALSE, missing=FALSE)
)
SPA_level_TF <- truef_spa %>%
select("Country", "BenchmarkID","PolicyActionImplemented") %>%
select(c("Country", "BenchmarkID")) %>% distinct() %>%
pivot_wider(names_from = "BenchmarkID", values_from = "BenchmarkID") %>%
ungroup() %>%
mutate(
across(.cols = -Country,
~if_else(!is.na(.), 1, 0)
)
)
AggHiearchy <- policydb %>% distinct(SubPArea, BenchmarkID, PolicyArea, PolicyAreaID)
PAreas <- policydb %>% distinct(PolicyArea, PolicyAreaID) #List of PAreas
SubPas <- policydb %>% distinct(SubPArea, BenchmarkID)
There are significant differences in number of policy actions collected per country:
In the following table it is evident that there are differences also by policy areas and sub-policy areas. The table is sort-able (e.g by PolicyArea).
So far, we know that there are differences in implemented of policy actions by both policy and sub-policy areas, as well as countries. This gives a clear indication that there will be differences in later stages of analyses.
The dataset is transformed to ensure consistency with the COINr framework. This is done to make a policy index where the policy dimensions are aggregated to a single number.
This an indicator of the indicators and units used the analyses. EEA countries are assigned to the EU, while EU countries not participating in the EU Comprehensive Scan outlined by WCRF is removed.
COINr_SPA <- SPA_level_TF
names(COINr_SPA)[names(COINr_SPA) == "Country"] <- "UnitName"
COINr_SPA <- COINr_SPA %>% mutate(UnitCode = countrycode(sourcevar = UnitName, origin="country.name", destination="iso3c"))
## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some values were not matched unambiguously: Gulf Cooperation Council, Micronesia, Romania, Wallis and Futuna
## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some strings were matched more than once, and therefore set to <NA> in the result: Romania, Wallis and Futuna,ROU,WLF
#Some values are missing: Gulf Cooperation Council, Micronesia, Romania, Wallis and Futuna
Missing <- COINr_SPA %>% subset(is.na(UnitCode)) %>% select(UnitName, UnitCode)
COINr_SPA <- na.omit(COINr_SPA)
isCoCREATE <- c("NOR", "NLD", "PRT", "POL", "GBR")
#Continue adding groups for countries
COINr_SPA <- COINr_SPA %>% mutate(
Group_Continent = countrycode(sourcevar = UnitName, origin="country.name", destination="continent"),
Group_EU = countrycode(sourcevar = UnitName, origin="country.name", destination="eu28"),
Group_Region = countrycode(sourcevar = UnitName, origin="country.name", destination="un.regionsub.name"),
Group_COCREATE = case_when(UnitCode %in% unlist(isCoCREATE) ~ TRUE, TRUE ~ FALSE))
## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some values were not matched unambiguously: Taiwan
COINr_SPA <- na.omit(COINr_SPA)
COINr_SPA <- COINr_SPA %>% mutate(Group_EU = case_when(
UnitCode == "NOR" ~ "EU",
UnitCode == "ISL" ~ "EU",
UnitCode == "LIE" ~ "EU",
UnitCode == "CHE" ~ "EU",
UnitCode == "RUS" ~ "NA", #not part of scan this round
TRUE ~ Group_EU)
)
COINr_SPA <- COINr_SPA %>% filter(Group_EU=="EU")
#only European countries for idnexation
COINr_SPA <- COINr_SPA %>% filter(Group_EU=="EU")
#Consider adding certain variables etc for further analysis (eurostat stuff)
IndData <- COINr_SPA
IndData <- IndData %>% mutate_if(is_double,as.integer)
IndData %>% reactable(compact = T)
The metadata for each indicators provides more details of single indicators and outlines the structure of the index.
IndMeta <- AggHiearchy
colnames(IndMeta) <- c("PolicyArea", "IndName", "IndCode", "Agg_PolicyDimension")
IndMeta$Direction <- 1 #All indicators got positive direction
IndMeta$IndWeight <- 1 #All indicators carry the same weight
#IndMeta$Agg1_Indicators <- IndMeta$IndCode ##Avventer denne. Blir dobbel uans?
#IndMeta$Agg_PolicyDimension #Samles til en av de tre hovedkategoriene
#IndMeta$Agg_PolicyIndex #Alt samles
IndMeta$Agg_PolicyIndex <- "NOURISHING"
IndMeta <- relocate(IndMeta, Agg_PolicyDimension, .before = Agg_PolicyIndex)
IndMeta %>% reactable(compact = T)
This outlines how the index aggregates from indicator to full index (3 steps).
AggMeta <- PAreas
colnames(AggMeta) <- c("Name", "Code")
AggMeta$AgLevel <- 2
AggMeta$Weight <- 1 #Equal weight for all
AggMeta <- AggMeta %>% add_row(Name="Index", Code="NOURISHING", AgLevel=3, Weight=1)
AggMeta %>% reactable(compact = T)
Combine the IndData, IndMeta and AggMeta into one hierarchical list. This creates a “COIN” and will be the base we work from here on out.
In this dataset we have 67 indicators (benchmarks), 31 units (countries) and two aggregation levels (indicators aggregate to a policy area and then to the full index).
NPI <- assemble(IndData = IndData, IndMeta = IndMeta, AggMeta = AggMeta)
NPI
## --------------
## A COIN with...
## --------------
## Input:
## Units: 32 (AUT, BEL, BGR, ...)
## Indicators: 67 (G1, G2, G3, ...)
## Denominators: 0 (none)
## Groups: 4 (Group_Continent, Group_EU, Group_Region, ...)
##
## Structure:
## Level 1: 67 indicators (G1, G2, G3, ...)
## Level 2: 10 groups (G, H, I, ...)
## Level 3: 1 groups (NOURISHING)
##
## Data sets:
## Raw (32 units)
The sunburst plot shows the hierarchical structure and effective weights for the simplified NOURISHING index.
Effective weights imply their relative contribution to the aggregated level, even with equal weighting applied. This is most apparent at indicator level: the N (level 2) has eight indicators contributing 0.0125 each, while in N(2), one indicator contributes 0.033. At sub-index level (N, O, U, R, I, S , H, I, N , G) they are contributing equally at 0.1 each.
An example of countries implementing R1:
The data indicate that some sub-policy areas are more covered than others (sort by “mean”).
The table indicate which areas have the ‘best’ coverage. This is I(2)1, N1 and N7 respectively. Multiple variables have low coverage (sort by MEAN), e. g. G5, H4 and R8.
Normalisation is not necessary if it is all on same scale, but moving from a smaller to larger scale will simplify presentatation of the results.
There are multiple ways to rescale data, but due to the binary nature of all indicators, a multiplication by 100 is sufficient. After this normalisation, the scores will either be 0 or 100 instead of 0-1.
NPI <- normalise(NPI, dset="Raw", ntype="custom", npara = list(custom = function(x) {x*100}))
A central step to create an index is the aggregation method. We do not want total compensability because being good in one area should not compensate for implementing less in another area.
In the first aggregation step, from indicator to policy area (e.g: from R1 to R) we use the arithmetic mean (simple averages).
The aggregation can be read left to right. The last column, NOURISHING, displays the total score for the index. Results are presented on the next page.
NPI <- aggregate(NPI, dset = "Normalised")
#A snapshot of the aggregation sequence
NPI$Data$Aggregated[(ncol(NPI$Data$Aggregated)-10): ncol(NPI$Data$Aggregated)] %>% roundDF(1) %>% reactable(compact = T)
The ranking of each country is in the plot and table below. Highlighted countries are members of the CO-CREATE project.
The countries perform differently according to policy areas as indicated by a deeper, green colour:
Examples using CO-CREATE countries. The overall mean for the index is (17.6475), but the implementation varies between the policy areas.
The spider charts are clickable. The first shows the overall mean, the second the mean of CO-CREATE countries.
iplotRadar(NPI, dset ="Aggregated", usel = isCoCREATE, aglev=2, addstat="mean")
iplotRadar(NPI, dset ="Aggregated", usel = isCoCREATE, aglev=2, statgroup = "Group_COCREATE", addstat = "groupmean")
The vast number of indicators and dimensions may make it hard to see patterns. The following analyses, we try to make two clusters: one for benchmarks and one for policy areas (as columns), with countries as the row.
#Prepare data
Aggregates <- getResults(NPI, tab_type = "Aggregates")
Aggregates <- Aggregates %>% mutate(
Group_Continent = countrycode(sourcevar = UnitName, origin="country.name", destination="continent"),
Group_EU = countrycode(sourcevar = UnitName, origin="country.name", destination="eu28"),
Group_Region = countrycode(sourcevar = UnitName, origin="country.name", destination="un.regionsub.name"),
Group_COCREATE = case_when(UnitCode %in% unlist(isCoCREATE) ~ TRUE, TRUE ~ FALSE))
Aggregates <- Aggregates %>% mutate(Group_EU = case_when(
UnitName == "NOR" ~ "EU",
UnitName == "ISL" ~ "EU",
UnitName == "LIE" ~ "EU",
UnitName == "CHE" ~ "EU",
TRUE ~ Group_EU)
)
#Data without text
NumericalsOnly <- Aggregates
rownames(NumericalsOnly) <- sapply(Aggregates$UnitName,function(x) strsplit(as.character(x),split = "\\\\")[[1]][1])
NumericalsOnly <- NumericalsOnly %>% select(-Rank, -NOURISHING, -UnitName, -UnitCode, -Group_EU, -Group_COCREATE, -Group_Region, -Group_Continent)
#Make unique group for visualisation
Regions <- Aggregates %>% select(Group_Region)
rownames(Regions) = rownames(NumericalsOnly)
The results for each policy area is set to divide into five categories of countries and three categories of policy areas.
NumericalsOnly2 <- NumericalsOnly %>% roundDF(decimals=1)
NumericalsOnly2[NumericalsOnly2 < 0.01] <- NA
NumericalsOnly2[is.na(NumericalsOnly2)] <- ""
library("RColorBrewer")
## Warning: package 'RColorBrewer' was built under R version 4.1.1
pheatmap::pheatmap(NumericalsOnly, cluster_cols = T, main="Clustering of countries by policy area", cutree_rows=3, cutree_cols=3, annotation_row = Regions, display_numbers = NumericalsOnly2, color=colorRampPalette(c("snow3", "yellowgreen", "forestgreen"))(50))
#Prepare data
AggregatesSubP <- getResults(NPI, tab_type = "Full")
#Data without text
NumericalsSuBPOnly <- AggregatesSubP
rownames(NumericalsSuBPOnly) <- sapply(NumericalsSuBPOnly$UnitName,function(x) strsplit(as.character(x),split = "\\\\")[[1]][1])
NumericalsSuBPOnly <- NumericalsSuBPOnly %>% select(-Rank, -NOURISHING, -UnitName, -UnitCode, -Group_EU, -Group_COCREATE, -Group_Region, -Group_Continent, -N, -O, -U, -R, -I, -S, -H, -"I.2.", -"N.2.", -G)
##ALTERNATIV
#Due to Pheatmap categories, I need to transpose the table to get the column annotation
trans_NumericalsSuBPOnly <- NumericalsSuBPOnly %>% t() %>% as.data.frame()
Dimensions <- AggHiearchy %>% select(BenchmarkID, PolicyAreaID) %>% as.data.frame()
Dimensions <- Dimensions %>% mutate(across("BenchmarkID", str_replace, "\\(", "."))
Dimensions <- Dimensions %>% mutate(across("BenchmarkID", str_replace, "\\)", "."))
#Add PolicyID and BenchmarkID
trans_NumericalsSuBPOnly <- left_join(trans_NumericalsSuBPOnly %>% mutate(BenchmarkID = rownames(trans_NumericalsSuBPOnly)), Dimensions, by = "BenchmarkID")
trans_NumericalsSuBPOnly <- trans_NumericalsSuBPOnly %>% column_to_rownames("BenchmarkID")
SubPAz <- trans_NumericalsSuBPOnly %>% select(PolicyAreaID)
NumericalsSuBPOnly <- trans_NumericalsSuBPOnly %>% select(-PolicyAreaID) %>% t() %>% as.data.frame()
When using the sub-policy areas to explore differences, the country ranking changes.
Important! These variables are “Yes” or “No”: hence the colour is either red for Implemented, or blue for Not implemented.
pheatmap::pheatmap(NumericalsSuBPOnly, cluster_cols = T, cluster_rows = T, main="Clustering of countries by benchmarks", cutree_rows=4, cutree_cols=3, annotation_col = SubPAz, legend = F, color = colorRampPalette(c("snow3", "forestgreen"))(50))
There are multiple approaches to using the final index score. An important prerequisite is a sound construction of the index. This will take some time to ensure and may end up not being feasible.
The regression can be used for a specific policy area (e. g., Nutrition labelling or Restricting marketing), or the full index. Both approaches should be explored. Note, there is less variance in some areas than others that we need to investigate closer.
We can retrieve any data from Eurostat or other sources. Here, I use Obesity statistics for the full population in European countries.
library(stringr)
hlth_ehis_bm1e <- eurostat::get_eurostat("hlth_ehis_bm1e")
## Table hlth_ehis_bm1e cached at C:\Users\JRMA\AppData\Local\Temp\RtmpMdZAhN/eurostat/hlth_ehis_bm1e_date_code_FF.rds
BMI_cat <- hlth_ehis_bm1e %>% select(bmi) %>% distinct()
hlth_ehis_bm1e$year <- str_sub(hlth_ehis_bm1e$time,1,4)
BMI_euro <- hlth_ehis_bm1e %>% filter(sex=="T" & age=="TOTAL" & isced11=="TOTAL" & bmi=="BMI_GE30" & year=="2019") %>% select(geo, values)
BMI_euro <- BMI_euro %>% mutate(geo2= countrycode(sourcevar = geo, origin="eurostat", destination="iso3c")) %>% select(geo2, values)
## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some values were not matched unambiguously: EU27_2020
colnames(BMI_euro) <- c("UnitCode", "Obesity rate")
BMI_euro %>% reactable::reactable()
After retriving the data from Eurostat, we combine the data from NOURISHING and the selected data source. Note, LIE has missing data. Data from different sources imported for GBR and CHE. For a final analysis, we should use HSBC or something else.
IndexScore_geo <- getResults(NPI, tab_type = "Summ") %>% select(UnitCode, NOURISHING)
IndexScore_geo <- left_join(IndexScore_geo, BMI_euro, by = c("UnitCode" ="UnitCode"))
#Add some missing data for this exercise. This data is not cross-checked!
IndexScore_geo <- IndexScore_geo %>% mutate("Obesity rate" = replace(`Obesity rate`, UnitCode=="GBR", 28))
IndexScore_geo <- IndexScore_geo %>% mutate("Obesity rate" = replace(`Obesity rate`, UnitCode=="CHE", 11.3))
IndexScore_geo %>% reactable::reactable(compact=T)
Creating a fitted line we find that there is a positive relationship between higher implementation of nutrition policies and obesity rates.
No significant correlation between index and obesity rates identified (p=0.35).
regline <- IndexScore_geo %>% filter(!is.na(`Obesity rate`)) %>% lm(NOURISHING ~ `Obesity rate`,.) %>% fitted.values()
corell <- cor.test(IndexScore_geo$NOURISHING, IndexScore_geo$`Obesity rate`, method = "pearson")
corell
##
## Pearson's product-moment correlation
##
## data: IndexScore_geo$NOURISHING and IndexScore_geo$`Obesity rate`
## t = 0.94849, df = 29, p-value = 0.3507
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1927261 0.4972377
## sample estimates:
## cor
## 0.1734597
x <- IndexScore_geo$NOURISHING
y <- IndexScore_geo$`Obesity rate`
Create a scatterplot with NOURISHING score on the vertical axis, and obesity rate on the horizontal axis:
meanNourish <- mean(IndexScore_geo$NOURISHING)
plot <- IndexScore_geo %>% filter(!is.na(`Obesity rate`)) %>%
plotly::plot_ly(x=~`Obesity rate`, y= ~NOURISHING, mode="markers", text =~UnitCode) %>%
plotly::add_markers(y= ~NOURISHING) %>%
plotly::add_trace(x=~`Obesity rate`, y=regline, mode="lines") %>%
plotly::add_lines(y=meanNourish) %>%
plotly::layout(showlegend=F) %>%
plotly::add_text(textposition="top right")
plot
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter